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The gravitational-wave (GW) sky may include nearby pointlike sources as well as astrophysical 
and cosmological stochastic backgrounds. Since the relative strength and angular distribution of the 
many possible sources of GWs are not well constrained, searches for GW signals must be performed 
in a model-independent way. To that end we perform two directional searches for persistent GWs 
using data from the LIGO S5 science run: one optimized for pointlike sources and one for arbitrary 
extended sources. The latter result is the first of its kind. Finding no evidence to support the 
detection of GWs, we present 90% confidence level (CL) upper-limit maps of GW strain power 
with typical values between 2 — 20 x 10 -50 strain 2 Hz _1 and 5 — 35 x 10 -49 strain 2 Hz -1 sr -1 for 
pointlike and extended sources respectively. The limits on pointlike sources constitute a factor of 
30 improvement over the previous best limits. We also set 90% CL limits on the narrow-band 
root-mean-square GW strain from interesting targets including Sco X-l, SN1987A and the Galactic 
Center as low as ~ 7 x 10 -25 in the most sensitive frequency range near 160 Hz. These limits are the 
most constraining to date and constitute a factor of 5 improvement over the previous best limits. 


I. INTRODUCTION 


One of the most ambitious goals of gravitational- wave 
(GW) astronomy is to measure the stochastic cosmo- 
logical gravitational- wave background (CGB), which can 
arise through a variety of mechanisms including ampli- 
fication of vacuum fluctuations following inflation [1], 
phase transitions in the early universe [2, 3], cosmic 
strings [4, 5] and pre-Big Bang models [6, 7]. The CGB 
may be masked by an astrophysical gravitational-wave 
background (AGB), interesting in its own right, which 
can arise from the superposition of unresolved sources 
such as core-collapse supernovae [8, 9], neutron-star ex- 
citations [10, 11], binary mergers [12, 13], persistent emis- 
sion from neutron stars [14, 15] and compact objects 
around supermassive black holes [16, 17]. 

We present the results of two analyses using data from 
the LIGO S5 science run: a radiometer analysis opti- 
mized for pointlike sources and a spherical-harmonic de- 
composition analysis, which allows for arbitrary angular 
distributions. This work presents the first measurement 
of the GW sky in a framework consistent with an arbi- 
trary extended source. 
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II. LIGO DETECTORS AND THE S5 SCIENCE 
RUN 

We analyze data from LIGO’s 4 km and 2 km detec- 
tors (HI and H2) in Hanford, WA and the 4 km detec- 
tor (LI) in Livingston Parish, LA during the S5 science 
run, which took place between Nov. 5, 2005 and Sep. 30, 
2007. During S5, both HI and LI reached a strain sen- 
sitivity of 3 x 10 -23 strain Hz -1 / 2 in the most sensitive 
region between 100 — 200 Hz [18] and collected 331 days 
of coincident H1L1 and H2L1 data. S5 saw milestones in 
GW astronomy including limits on the emission of GWs 
from the Crab pulsar that surpass those inferred from the 
Crab’s spindown [19], as well as limits on the isotropic 
CGB that surpass indirect limits from Big Bang nucle- 
osynthesis and the cosmic microwave background [20]. 
This work builds on [20, 21]. 


III. METHODOLOGY 

Following [21, 22] we present a framework for analyz- 
ing the angular distribution of GWs. We assume that the 
GW signal is stationary and unpolarized, but not neces- 
sarily isotropic. It follows that the GW energy density 
9 gw (/), can be expressed in terms of the GW power 
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spectrum, V(f,Cl): 

^GW (/) = — = 17772 f 3 [ dClV(f,Cl). (1) 

p c aj 6H 0 J s 2 

Here / is frequency, Cl is sky location, p c is the criti- 
cal density of the universe and Hq is Hubble’s constant. 
We further assume that V(f,Cl) can be factored (in our 
analysis band) into an angular power spectrum, V(Cl ), 
and a spectral shape, H(f ) = (///o)^, parameterized by 
the spectral index (3 and reference frequency /o- We set 
/o = 100 Hz to be in the sensitive range of the LIGO 
interferometers . 

Our goal is to measure V{Cl) for two power-law signal 
models. In the cosmological model, [3 = —3 (Dqw(/) 
const), which is predicted, e.g., for the amplification of 
vacuum fluctuations following inflation (see [23] and ref- 
erences therein). In the astrophysical model, (3 = 0 
(H(f) = const), which emphasizes the strain sensitivity 
of the LIGO detectors. 

We estimate V(Cl) two ways. The radiometer algo- 
rithm [21, 24, 25] assumes the signal is a point source 
characterized by a single direction flo and amplitude, 
rj(£l 0 ): 

v{h) = v {n 0 )s 2 (Ci,Ci 0 ). ( 2 ) 

It is applicable to a GW sky dominated by a limited 
number of widely separated point sources. As the number 
of point sources is increased, however, the interferometer 
beam pattern will cause the signals to interfere and partly 
cancel. Thus, radiometer maps do not apply to extended 
sources. Since pointlike signals are expected to arise from 
astrophysical sources, we use j3 = 0 for the radiometer 
analysis. 

The spherical-harmonic decomposition (SHD) algo- 
rithm is used for both (3 = — 3 (cosmological) and /3 = 0 
(astrophysical) sources. It allows for the possibility of an 
extended source with an arbitrary angular distribution, 
characterized by spherical-harmonic coefficients Vi m such 
that 

r(n) = J2'Pi m Y lm ^)- ( 3 ) 

Im 

The series is cut off at / max , allowing for angular scale 
~ 2tt // max . The flexibility of the spherical-harmonic al- 
gorithm comes at the price of somewhat diminished sen- 
sitivity to point sources, and thus the two algorithms are 
complementary. 

We choose / max so as to minimize the sky average of 
the product of cr(D) and A, where cr(II) is the uncertainty 
associated with V(Cl) and A is the typical angular area 
of a resolved patch of sky [37]. Since A = 47r/7Vi n d ep oc 
l/(GaxT l) 2 (where Wndep is the number of independent 
parameters), this procedure amounts to choosing / max to 
maximize the sensitivity obtained by integrating over the 
typical search aperture (angular resolution). We obtain 
imax = 7 and 12 for /3 = —3 and (3 = 0, respectively. 


Since the search aperture becomes smaller at the higher 
frequencies emphasized by (3 = 0 , / max is larger for (3 = 0 
than for f3 = —3. 

Both algorithms can be framed in terms of a “dirty 
map”, X u , which represents the signal convolved the 
Fisher matrix, T Ml/ : 

X " = E 7; i (/, t) Pi (/ , t) C(f, t) (4) 

Here both the Greek indices fi and v take on values of 
Im for the SHD algorithm and Cl for the radiometer al- 
gorithm, for which we use the pixel basis. The two bases 
are related using spherical-harmonic basis functions: 

X Cl = J2 X lmYlm^)- ( 6 ) 

Im 

C(/, £), meanwhile, is the cross spectral density gener- 
ated from the H1L1 or H2L1 pairs. Pi(/, t) and P 2 (f,t) 
are the individual power spectral densities, and 7 M (/, t ) is 
the angular decomposition of the overlap reduction func- 
tion 7 (D, /,t), which characterizes the orientations and 
frequency response of the detectors [ 22 ]: 

74/^ ) = [ d£l'y(£l,f,t)e ll (£l) (7) 

Js 2 

7 = ^F^(n,t)F^(n,t)e i27TfCl < ASl2 ^ t)) / c . ( 8 ) 

Here Ff(Q,t) characterizes the detector response of de- 
tector / to a GW with polarization A , e^(Cl) is a basis 
function, c is the speed of light and Axi 2 = x\ — X 2 is 
the difference between the interferometer locations. A de- 
tailed discussion of these quantities can be found in [ 22 ]. 

In [22] it was shown that the maximum-likelihood es- 
timators of GW power are given by V — Y~ 1 X. The 
inversion of T is complicated by singular eigenvalues asso- 
ciated with modes to which the Hanford-Livingston (HL) 
detector network is insensitive. This singularity can be 
handled two ways. The radiometer algorithm assumes 
the signal is pointlike, implying that correlations between 
neighboring pixels can be ignored. Consequently, we can 
replace T _1 with (T ^) -1 to estimate the point source 
amplitude r](Cl) (see Eq. 2 ). (We note that pointlike 
sources create signatures in our sky maps that typically 
span several degrees or more; see [ 21 ].) 

The SHD algorithm, on the other hand, targets ex- 
tended sources, so the full Fisher matrix must be taken 
into account. We regularize T by removing a fraction, 
P, of the modes associated with the smallest eigenval- 
ues, to which the HL network is relatively insensitive. T 
is known as the regularization cutoff. By removing some 
modes from the Fisher matrix, we obtain a regularized 
inverse Fisher matrix, T^ 1 , thereby introducing a bias, 
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the implications of which are discussed below. For now, 
we note that the bias depends on the angular distribution 
of the signal. 

We thereby obtain the estimators 

% = ft tin) (9) 

Vim = ^(r R 1 )lm,l'm'Xi'm' , (10) 


with uncertainties 


bins up to 200,500 Hz, so that a (Cl) is within < 2% of 
the minimum possible value. Thirty-three frequency bins 
are masked, corresponding to 2% of the frequency bins 
between 40 — 500 Hz used in the broadband analyses. 
For additional details about the S5 stochastic pipeline, 
see [20]. 


IV. SIGNIFICANCE AND UPPER LIMIT 
CALCULATIONS 


^ ad = (ran)" 172 (11) 

a lm = [(r R)lm,lm\ • (12) 

We refer to Vq = Y^im^imYim(Cl) as the “clean map” 
and fjfi as the “radiometer map.” We note that 
fj ^ has units of strain 2 Hz -1 whereas has units of 
strain 2 Hz -1 sr -1 . 

In choosing T one must balance the competing de- 
mands of reconstruction accuracy (sensitivity to the 
modes that are kept) with the bias associated with the 
modes that are removed. In practice, we do not know the 
bias associated with T since it depends on the unknown 
signal distribution V(Cl). Therefore, we choose a value of 
T that tends to produce reliably reconstructed maps with 
minimal bias for simulated signals. Following [22], we use 
T — 1/3, which was shown to be a robust regularization 
cutoff for simulated signals including maps characterized 
by one or more point sources, dipoles, monopoles and an 
extended source clustered in the galactic plane (see [22]). 

In the case of the SHD algorithm, we construct an 
additional statistic (see [22]), 


Ci = 


2Z + 1 ^ 


[l A 


(r 




(13) 


which describes the angular scale of the clean map. The 
subtracted second term makes the estimator unbiased so 
that (Ci) = 0 when no signal is present. The expected 
noise distribution of Ci is highly non- Gaussian for small 
values of /, and so the upper limits presented below are 
calculated numerically. The Ci are analogous to similar 
quantities defined in the context of temperature fluctua- 
tions of the cosmic microwave background (see, e.g., [26]). 

The analysis was performed using the S5 stochastic 
analysis pipeline. This pipeline has been tested with 
hardware and software injections, and the successful re- 
covery of isotropic hardware injections is documented 
in [20]. The recovery of anisotropic software injections 
is demonstrated in [22]. We parse time series into 60s, 
Hann- windowed, 50%-overlapping segments, which are 
coarse-grained to achieve 0.25 Hz resolution. We ap- 
ply a stationarity cut described in [21], which rejects 
^ 3% of the cross-correlated segments. We also mask 
frequency bins associated with instrumental lines (e.g., 
harmonics of the 60 Hz mains power, calibration lines 
and suspension-wire resonances) as well as injected, sim- 
ulated pulsar signals. For {3 = —3, 0 we include frequency 


In order to determine if there is a statistically signif- 
icant GW signature, we are primarily interested in the 
significance of outliers — the highest signal-to-noise ratio 
(SNR) frequency bin or sky-map pixel. It is therefore 
necessary to calculate the expected noise probability dis- 
tribution of the maximum SNR given many independent 
trials (when considering maximum SNR in a spectral 
band) and given many dependent trials (when considering 
maximum SNR for a sky map). 

For N independent frequency bins, the probability 
density function, 7r(p max ), of maximum SNR, p max , is 
given by 


7r(/?max) OC 


1 + erf p n 


N-l 


(14) 


Here we have assumed that the stochastic point esti- 
mate is Gaussian distributed. The Gaussianity of and 
t)^, calculated by summing over many 0(5OOK) indepen- 
dent segments, is expected to arise due to the central 
limit theorem [27]. Additionally, we find the Gaussian- 
noise hypothesis to be consistent with time-slide studies, 
wherein we perform the cross-correlation analysis with 
an unphysical time-shift in order to wash out astrophys- 
ical signals and thereby obtain different realizations of 
detector noise. 

The distribution of maximum SNR for a sky map is 
more subtle due to the non-zero covariances that exist 
between different pixels (or patches) on the sky. For 
this case, we calculate 7r(p max ) numerically, by simu- 
lating many realizations of dirty maps that have ex- 
pected covariances described by the Fisher matrix T. Fig- 
ure 1 shows the numerically determined 7r(p max ) for the 
/3 = —3 clean map generated with Gaussian noise. 

The likelihood function for V(Cl) at each point in the 
sky can be be described as a normal distribution with 
mean 'Pn and a variance (cr^ ph ) 2 . In the case of the SHD 
algorithm, regularization introduces a signal-dependent 
bias. Without knowing the true distribution of 7^(0), it 
is impossible to know the bias exactly, but it is possible 
to set a conservative upper limit by assuming that on av- 
erage the modes removed through regularization contain 
no more GW power than the modes that are kept. 

To implement this assumption, we calculate Vi m with 
a regularization scheme that sets eigenvalues of removed 
modes to zero, whereas is conservatively calculated 
using a regularization scheme that sets eigenvalues of 
removed modes to the average eigenvalue of the kept 
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FIG. 1: Numerically calculated distribution of the maximum 
SNR for f3 = — 3 clean maps created from Gaussian noise. 


modes. This has the effect of widening the likelihood 
function at each sky location. The upper limits become 
on average 25% larger than they would be if we had calcu- 
lated cr^ h using the same regularization scheme as Vi m . 

Following the same procedure as in [20], we marginal- 
ize over the HI, H2, and LI calibration uncertainties, 
which were measured to be 10%, 10%, and 13%, re- 
spectively [28] [38] The posterior distribution is obtained 
by multiplying the marginalized likelihood function by a 
prior, which we take to be flat above V{Ct) > 0 [39]. The 
Bayesian upper limits are then determined by integrat- 
ing the posterior out to the value of V(Cl) which includes 
90% of the total area under the distribution. The calcu- 
lation of upper limits on 77^ is analogous except we need 
not take into account the effects of regularization. 


V. RESULTS 


Sky maps: Figure 2 presents sky map results for the 
different analyses: SHD algorithm with (3 = — 3 (left), 
SHD with (3 = 0 (center), and radiometer with (3 = 0 
(right). The top row contains SNR maps. The maximum 
SNR values are 3.1 (with significance p = 25%), 3.1 (with 
p = 56%), and 3.2 (with p = 53%) respectively. These 79- 
values take into account the number of search directions 
and covariances between different sky patches (see IV). 
Observing no evidence of GWs, we set upper limits on 
GW power as a function of direction. The 90% confidence 
level (CL) upper limit maps are given in the bottom row. 
For the SHD method with f3 = —3, the limits are between 
5 — 31 x 10 -49 strain 2 Hz _1 sr _1 ; for SHD with (3 = 0, the 
limits are between 6 — 35 x 10 -49 strain 2 Hz -1 sr -1 ; and 
for the radiometer with (3 = 0, the limits are between 
2-20 x 10 -5 ° strain 2 Hz _1 . 

The strain power limits can also be expressed in terms 


of the GW energy flux per unit frequency [21]: 

= ( 3 18 x 1042 S) (iss ) 9+2 < 16 > 

(Radiometer energy flux is obtained by replacing 
with fjft.) The corresponding values are 2 — 

10 x 10 _6 (//100 Hz) -1 ergcm _2 s _1 Hz _1 sr _1 and 2 — 

11 x 10 _6 (//100 Hz) 2 ergcm _2 s _1 Hz _1 sr _1 for the SHD 
method, and 6 — 60 x 10 _8 (//100 Hz) 2 ergcm _2 s _1 Hz _1 
for the radiometer. The radiometer map constitutes a 
factor of ~ 30 improvement over the previous best strain 
power limits [21]. 

When comparing the SHD analysis with /3 = 0 and the 
radiometer upper limits obtained using the same spec- 
trum, it is important to note that these maps have differ- 
ent units. The radiometer map has units of strain 2 Hz -1 
because the radiometer analysis effectively integrates the 
power from a GW point source over solid angle. The SHD 
maps, on the other hand, have units of strain 2 Hz _1 sr _1 . 
If we scale the SHD limit maps by the typical diffraction 
limited resolution (A ~ O.lsr), then the limits are more 
comparable. The radiometer algorithm limits are lower 
(by a factor of < 2) because it requires a stronger as- 
sumption about the signal model (a single point source), 
whereas the SHD algorithm is model- independent. 

Figure 3 show 90% CL upper limits on the Ci. Since 
the Vim Lave units of strain power (strain 2 Hz -1 sr -1 ), the 
Ci have the somewhat unusual units of strain 4 Hz -2 sr -2 . 

Targeted searches: Sco X-l is a nearby (2.8 kpc) low- 
mass X-ray binary likely to include a neutron star spun 
up through accretion. Its spin frequency is unknown. It 
has been suggested that this accretion torque is balanced 
by GW emission [29]. The Doppler line broadening due 
to the orbital motion is smaller than the chosen Sf = 
0.25 Hz bin width for frequencies below ~ 930 Hz [30]. 
At higher frequencies, the signal is certain to span two 
bins. We determine the maximum value of SNR in the di- 
rection of Sco X-l to be 3.6 at / = 1770.50 Hz, which has 
a significance of p = 73% given the 0(7000) independent 
frequency bins. Thus in Fig. 4 (first panel) we present 
limits on root-mean-square (RMS) strain, /irms(/>^)> 
as a function of frequency in the direction of Sco X-l 
(RA, dec) = (16.3 hr, 15.6°). These limits improve on the 
previous best limits by a factor of ^ 5 [21]. RMS strain 
is related to narrow-band GW power via 


hRMs(fP)= r)(f,£l)Sf 


1 1/2 


(17) 


and is better suited for comparison with searches for pe- 
riodic GWs, which typically constrain the peak strain 
amplitude, Lo, marginalized over neutron star parame- 
ters l and ^ (see, e.g., [31]). Our limits on Lrms are 
for a circularly polarized signal from a pulsar whose spin 
axis is aligned with the line of sight. Marginalizing over 
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FIG. 2: Top row: Signal-to-noise ratio maps for the three different analyses described in this paper: SHD clean map for f3 = — 3 
(left), SHD clean map for ft = 0 (center), and radiometer for /3 = 0 (right). All three SNR maps are consistent with detector 
noise. The p-values associated with each map’s maximum SNR are (from left to right) p = 25%, p = 56%, p m 53%. Bottom 
row: The corresponding 90% CL upper limit maps on strain power in units of strain 2 Hz - 1 sr _1 for the SHD algorithm, and 
units of strain 2 Hz -1 for the radiometer algorithm. 
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FIG. 3: Upper limits on Ci at 90% CL vs l for the SHD analyses for [3 = — 3 (left) and /3 = 0 (right). The Ci are consistent 
with detector noise. 


l and ip and converting from RMS to peak amplitude 
causes the limits to change by a sky-dependent factor of 
~ 2.3 [32]. We note that these limits are on the RMS 
strain in each bin as opposed to the total RMS strain 
from Sco X-l, which might span as many as two bins. 
The frequency axis refers to the observed GW frequency 
as opposed to the intrinsic GW frequency. 

We also look for statistically significant outliers 
associated with the Galactic Center (RA, dec) = 
(17.8 hr, -29°) and SN1987A (RA,dec) = (5.6 hr, -69°). 
The maximum SNR values are 3.5 at / = 203.25 Hz with 
p = 85% and 4.3 at 1367.25 Hz with p = 7%, respec- 
tively. Limits on RMS strain are given in the right panel 
of Fig. 4. 


VI. CONCLUSIONS 

We performed two directional analyses for persis- 
tent GWs: the radiometer analysis, which is optimized 
for point sources, and the complementary spherical- 
harmonic decomposition (SHD) algorithm, which allows 
for arbitrary extended sources. Neither analysis finds evi- 
dence of GWs. Thus we present upper-limit maps of GW 
power and also limits on the RMS strain from Sco X-l, 
the Galactic Center and SN1987A. The radiometer map 
limits improve on the previous best limits [21] by a fac- 


tor of 30 in strain power, and limits on RMS strain from 
Sco X-l constitute a factor of 5 improvement in strain 
over the previous best limits [21]. The SHD clean maps 
represent the first effort to look for anisotropic extended 
sources of GWs. 

With the ongoing construction of second-generation 
GW interferometers, we are poised to enter a new era 
in GW astronomy. Advanced detectors [33-36] are ex- 
pected to achieve strain sensitivities approximately 10 
times lower than initial LIGO, and advances in seismic 
isolation are expected to extend the frequency band down 
from 40 Hz to 10 Hz [33]. By adding additional detec- 
tors to our network, we expect to reduce degeneracies in 
the Fisher matrix and improve angular resolution. These 
improvements will allow advanced detector networks to 
probe plausible models of astrophysical stochastic fore- 
grounds and some cosmological models such as cosmic 
strings. 
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